Simulating Impacts of Precipitation on Ice Cover and Surface Water Temperature Across Large Lakes
Abstract
Precipitation impacts on ice cover and water temperature in the Laurentian Great Lakes were examined using state-of-the-art coupled ice-hydrodynamic models. Numerical experiments were conducted for the recent anomalously cold (2014–2015) and warm (2015–2016) winters that were accompanied by high and low ice coverage over the lakes, respectively. The results of numerical experiments showed that snow cover on the ice, which is the manifestation of winter precipitation, reduced the total ice volume (or mean ice thickness) in all of the Great Lakes, shortened the ice duration, and allowed earlier warming of water surface. The reduced ice volume was due to the thermal insulation of snow cover. The surface albedo was also increased by snow cover, but its impact on the delay the melting of ice was overcome by the thermal insulation effect. During major snowstorms, snowfall over the open lake caused notable cooling of the water surface due to latent heat absorption. Overall, the sensible heat flux from rain in spring and summer was found to have negligible impacts on the water surface temperature. Although uncertainties remain in overlake precipitation estimates and model's representation of snow on the ice, this study demonstrated that winter precipitation, particularly snowfall on the ice and water surfaces, is an important contributing factor in Great Lakes ice production and thermal conditions from late fall to spring.
Key Points
- Precipitation impacts on Great Lakes ice cover and water temperature were evaluated using a coupled ice-hydrodynamic model
- The model results showed that snow cover on the ice reduced the net production of ice, which resulted in earlier decay of ice cover
- The model results showed that snowfall cooled the water surface notably through latent heat absorption during storms
Plain Language Summary
Snow and rain impact on ice cover and water temperature in large lakes were studied using a computational model for an example of the Laurentian Great Lakes. It was found that snow cover increased the reflection of solar radiation but at the same time prevented lake ice from the growing, resulting in less formation of ice and slightly earlier melting. The earlier ice melting also allowed earlier warming of the water surface in spring. Major snowstorms caused slight cooling in the water surface temperature because snowflakes absorbed heat when it touched the water surface to melt. On the other hand, warmer rain barely changed the water surface temperature during summer.
1 Introduction
Understanding precipitation dynamics over Earth's large freshwater surfaces is critical in order to reconcile fluxes of energy and moisture across global and continental scales (National Research Council, 2007). While the role of precipitation has been recognized as a major factor in water balance from a hydrology perspective, precipitation impacts on ice and water temperature across Earth's large lakes are relatively undocumented. Understanding these relationships becomes all-the-more important as physical properties of lakes evolve in response to climate change and water management practices (Erler et al., 2019; Khazaei et al., 2019; Wuebbles et al., 2019); in some cases, these changes have led to severe lake water level declines and even lake disappearance (Gao et al., 2011; Shi et al., 2019; Tao et al., 2015).
In middle- and high-latitude lakes, there are a few rationales for why precipitation can be important in these processes. First, snow accumulation on lake ice, which is a manifestation of winter precipitation, has two opposing effects on lake ice. It increases the ice surface albedo, reflecting more incoming shortwave radiation than bare ice (Gardner & Sharp, 2010; Wiscombe & Warren, 1980), and slows down ice melting from the surface. On the other hand, because the heat conductivity of snow is generally lower than that of lake ice (Quinn et al., 1978; Yen, 1981), it insulates ice underneath from atmospheric cooling and therefore slows down the ice growth from the ice-water interface. These two opposing impacts are well known (Fichefet & Morales Maqueda, 1999; Sturm et al., 2002; Warren et al., 1999), but empirical studies and model simulations to assess their impacts on sea ice or lake ice are rare, partly because of the high uncertainty in overwater precipitation (Holman et al., 2012) and the complexity of snow morphology (Webster et al., 2018).
Second, the air-lake heat transfer associated with precipitation can be significant. This heat transfer can be divided into two components, that is, the sensible and latent components. The sensible heat flux from precipitation occurs due to the temperature difference between rain droplets/snowflakes and the lake surface. There are a limited number of studies on this heat transfer in oceans (e.g., Duffy & Bennartz, 2018; Gosnell et al., 1995) and over glaciers (e.g., Alexander et al., 2011; Anderson et al., 2010), most of which concluded that the sensible heat flux from precipitation is negligible compared with the other major components of the heat budget. However, these studies focused on places where the temperature gradient between the atmosphere above and the surface was relatively low (<5°C). However, the situation is different in middle- and high-latitude lakes. In the Laurentian Great Lakes (hereafter Great Lakes); for example, the atmosphere-lake temperature can be large (>10°C), especially during the fall, and therefore, the sensible heat flux from precipitation is larger. The latent heat flux from precipitation occurs due to melting of snow when it touches down on the lake surface. Unlike the tropical ocean, massive snowstorms over the Great Lakes and any large lake in middle and high latitudes during winter may cause significant latent heat flux due to snow melting. This process is gaining attention and has been implemented in numerical models for high-latitude oceans (e.g., Duffy & Bennartz, 2018; Durski & Kurapov, 2019); thus, a curiosity exists regarding the importance of snow melting in the Great Lakes heat budget.
There is a growing momentum in the coastal modeling community for coupling ice, hydrodynamics, and hydrologic processes (Elko et al., 2014; Liu et al., 2019; Munar et al., 2018). Examining precipitation impacts on the Great Lake ice and water temperature would be a suitable contribution to ensuring accurate interactions at the lake surface in coupled model applications. In this study, we examine the impacts of precipitation on the seasonal evolution of lake ice and water temperature in the Great Lakes, a prime example of midlatitude large lakes, using a state-of-the-art ice-hydrodynamic model, driven by a forcing data set from a high-resolution weather forecast model.
In section 2, we describe the model and data sets used in the study. In section 3, we present the results from the model simulations as well as comparison with available observations. In section 3, also, we discuss the implications of the results and the importance of precipitation impacts to seasonal ice simulation, as well as on air-lake coupling in weather forecast applications. In section 4, we summarize and conclude our findings.
2 Materials and Methods
2.1 Ice-Hydrodynamic Model
We use the unstructured grid, Finite Volume Community Ocean Model (FVCOM, Chen et al., 2006, 2013) to simulate the Great Lakes hydrodynamics. FVCOM is a three-dimensional, free-surface, primitive equation, sigma-coordinate oceanographic model that solves the integral form of the governing equations. FVCOM has been applied in several studies of the coastal ocean, including successful application to operational forecasting in the Great Lakes (E. J. Anderson et al., 2010, 2015, 2018; E. J. Anderson & Schwab, 2012, 2013; Bai et al., 2013; Niu et al., 2015; Luo et al., 2012). In this work, the model is configured separately for Lake Superior, Lake Erie, and Lake Ontario, while Lake Michigan and Lake Huron are handled by the same model as the lakes are connected by the Straits of Mackinac and form a single system. Horizontal grid resolution in each model ranges from roughly 200 m near the shoreline to 2,500 m offshore, with 21 vertical sigma layers evenly distributed throughout the water column. As a result, the numbers of triangular elements in the models are roughly 23,000 for Lake Superior, 31,000 for Lake Michigan-Huron, 12,000 for Lake Erie, and 35,000 for Lake Ontario. For Lake Erie, the grid configuration is based on the National Oceanic and Atmospheric Administration (NOAA)'s Lake Erie Operational Forecast System (Anderson et al., 2018), and the remaining lakes have similar resolution and configuration. Horizontal diffusion is handled by the Smagorinsky parameterization (Smagorinsky, 1963), and vertical diffusion is handled by Mellor-Yamada Level-2.5 turbulence closure scheme (Mellor & Blumberg, 2004; Mellor & Yamada, 1982) with surface wave breaking parameterization by Craig and Banner (1994). Turbulent latent and sensible heat fluxes (not associated with precipitation) are calculated from the Coupled Ocean-Atmosphere Response Experiment (COARE, Fairall et al., 1996, Fairall et al., 1996, Fairall et al., 2003) algorithm. Modeled depths (Figure 1) are interpolated from 3-arc sec bathymetry data from the NOAA National Centers for Environmental Information.

An unstructured grid version of the Los Alamos Sea Ice model (CICE; Gao et al., 2011; Hunke et al., 2015) has been coupled within FVCOM. The CICE model includes components for ice thermodynamics and ice dynamics, using elastic-viscous-plastic rheology for internal stress (E. C. Hunke & Dukowicz, 1997) and produces two-dimensional fields of ice concentration, thickness, and velocity. A multicategory ice thickness distribution model (Thorndike et al., 1975) is employed in CICE to represent the sub-grid-scale distribution of ice thickness in response to mechanical and thermal forcing. In this study, five categories of ice thickness are defined (5, 25, 65, 125, and 205 cm). The ice surface albedo depends on surface temperature, snow depth, and thickness of ice, as well as the visible and infrared spectral bands of the incoming solar radiation (Briegleb, 1992). Specifically, visible and near-infrared spectral albedos are distinguished for snow and ice. For the bare ice (snow) albedo, 0.78 (0.98) and 0.36 (0.70) are used for visible (<700 nm) and near-infrared (>700 nm) wavelengths, respectively. The ice (snow) albedos in both spectral bands decrease by 0.075 (0.15) as the surface temperature rises from −1°C to 0°C. As ice thickness decreases from 50 cm to 0, the ice albedo decreases smoothly (using an arctangent function) to the open water albedo, 0.06. The total albedo is an area-weighted average of the ice and snow albedos, where the fractional snow-covered area fsnow is defined as fsnow = hs/(hs + hsnowpatch), where hs is snow depth and hsnowpatch is a constant set at 2 cm. The snow and ice albedo parameterization is consistent with the one in the previous version of the Community Climate System Model (CCSM3) and is detailed in Hunke et al. (2019). For the density of snow ρs, a constant value of 330 kg/m3 was used. The thermal conductivity of snow and ice were set to 2.03 and 0.30 W/m/K, respectively.
At ice-covered cells, the net momentum transfer is calculated as a weighted average of the air-water and ice-water stresses by areal fraction of ice. The air-ice drag coefficient CD_ai is a function of wind speed U, given as CD_ai = (1.43 + 0.052 U) × 10−3 and the ice-water drag coefficient is 5.5 × 10−3. Similarly, the net heat transfer is calculated as a weighted average of the air-water and ice-water heat fluxes. The ice-water heat fluxes are calculated based on the bulk transfer formula (Maykut & McPhee, 1995). The formation of snow ice, which forms when flooded snow refreezes, is parameterized by calculating the depth of snow submerged below the water surface based on the hydrostatic assumption (CICE; Consortium, 2019). The submerged snow is converted to ice.
The selected experiment period was from 11 November 2014 to 1 November 2016. This period covers one severe winter with high ice cover (2014–2015) and a mild winter with low ice cover (2015–2016) and is ideally suited to evaluate the impacts of precipitation on variable ice conditions, as well as thermal conditions afterwards. The control case started on 1 April 2014 with uniform 4°C temperature and static conditions (cold start option) and forced by interpolated observations from coastal and offshore meteorological stations that are corrected for overwater conditions (Beletsky et al., 2003) until September 2014, after which an operational weather forecast from the High-Resolution Rapid Refresh (HRRR, described in section 2.2) was utilized as the forcing data. All the three experiments, which are described in section 2.3, were started from a restart file of this control case on 11 November 2014. This start date for the three numerical experiments allowed to cover the first winter storm over the Great Lakes region in the winter of 2014–2015 (Fujisaki-Manome et al., 2017) when notable snow and therefore was appropriate to examine the impacts of precipitation (snow or rain) on lake ice and water temperature.
2.2 Data
For the meteorological forcing to the FVCOM-CICE model, the HRRR (Benjamin et al., 2016,Benjamin et al., 2016) was selected. HRRR is a relatively new, cloud-resolving, convection-allowing weather model with horizontal resolution of 3 km that became operational in September 2014. Radar data are assimilated in the HRRR every 15 min over a 1-hr period, adding further details to those provided by the hourly data assimilation from the 13-km radar-enhanced Rapid Refresh. The meteorological forcings from HRRR at forecast hour 2 were applied to the FVCOM-CICE model on an hourly basis by providing wind speed at 7.8-m height, surface air temperature at 2-m height, relative humidity at 2-m height, and precipitation rate. The COARE algorithm used in this study is able to take account of these reference heights to derive the momentum and turbulent heat fluxes in the surface boundary layer.
Cross-evaluation of precipitation rate from the HRRR was conducted beforehand in comparison with the merged product of the Canadian Precipitation Analysis (CaPA) and the Multisensor Precipitation Estimate (MPE). CaPA is a real-time gridded precipitation product provided by Environment and Climate Change Canada and covers all of North America. MPE is the composite analysis of weather radars, rain gage data, and satellite precipitation estimates and is provided by the U.S. National Weather Service. The merged CaPA-MPE product combines these two data sets into 10-km gridded data and has spatially continuous precipitation estimate over the Great Lakes. The merged CaPA-MPE product is widely used by the binational water balance community for hydrologic analysis and lake management and therefore an adequate data set to compare with HRRR's overlake precipitation to ensure its validity. Further details can be referred to Gronewold et al. (2018). We calculated daily overlake precipitation for each lake from the HRRR and compared with those from the merged CaPA-MPE product (Figure 2). The two data sets generally agreed well. Slightly larger values in HRRR, represented by a slope value s > 1, are likely because the CaPA-MPE rely heavily on radar estimates from the coastlines, which are sometimes limited in terms of coverage and obstacles (e.g., the radar at Marquette, Michigan, subject to the blockage by the Keweenaw Peninsula), and therefore could miss some overwater precipitation. Overall, the reasonable agreement of HRRR's overlake precipitation estimate with the merged CaPA-MPE provides confidence in using overlake precipitation from HRRR. It should be reminded that precipitation from HRRR has other advantages such as high spatiotemporal resolution and physically consistent representation based on the single model.

For the evaluation of water surface temperature and ice coverage, the Great Lakes Surface Environmental Analysis (GLSEA; https://coastwatch.glerl.noaa.gov/glsea/doc/) was used along with ice concentration maps from the National Ice Center. GLSEA is the composite analysis from six satellite measurements based on the Advanced Very High Resolution Radiometer and the Visible Infrared Imaging Radiometer Suite. The water surface temperatures are updated daily using information from the cloud-free portions of the previous day's satellite imagery. For pixels where data are not available, an interpolation algorithm is applied using cloud-free data over a ±10-day window. Therefore, the accuracy tends to be lower on cloudy days.
2.3 Numerical Experiments
We conducted three numerical experiments to assess the impacts of precipitation on lake ice and water temperature. In the control experiment (Expt. 1), no precipitation is considered. In the first precipitation experiment (Expt. 2), precipitation rate from the HRRR was applied to the model as the input data. The sensible and latent heat fluxes from precipitation (Hsp and Hlp) were calculated by the models in these experiments. The other heat flux components (i.e., shortwave radiation, longwave radiation, turbulent sensible, and latent heat fluxes) were calculated dynamically in the model and Hsp and Hlp are aggregated with the other heat flux components into the net heat flux. We also conducted the supplemental precipitation experiment (Expt. 3), where the other heat flux components than Hsp and Hlp were prescribed from the outputs of Expt. 1. The net heat flux over water was extracted from the Exp. 1 (no precipitation) results, and this net heat flux was fed to Expt. 3, where the model calculated the heat fluxes only for Hsp and Hlp and added them to the prescribed net heat flux. The purpose of Expt. 3 was to evaluate the maximum potential of Hsp and Hlp impacts on the water temperature. In Expt. 2, both water temperature and the other heat flux components were allowed to respond to the precipitation heat fluxes (i.e., perturbation), while in Expt. 3, only water temperature was allowed to respond to the perturbation because the other heat flux components were fixed. Therefore, the impacts on water temperature are meant to be maximized in Expt. 3.
3 Results and Discussions
3.1 Verification of Modeled Ice Conditions
Seasonal evolutions of ice coverage were reasonably simulated for each of the Great Lakes by the models (Figure 3) in comparison with the National Ice Center (NIC) analysis. This is consistent with the long-term model validation for Lake Erie and Lake Michigan-Huron by Anderson et al. (2018). The contrasts of high and low maximum ice coverage reflected the anomalously cold (2014–2015) and warm (2015–2016) winters. The spatial patterns of ice concentration at the early, middle, and late ice stages agreed with the NIC analysis (Figure 4). There were some overestimations of ice coverage during late ice season (March–April) in 2015, which were most notable in Lakes Superior and Huron and in Lakes Erie and Ontario (Figure 3). With precipitation (Expts. 2 and 3), ice coverage declined slightly faster than that in Expt. 1. The root-mean-square errors (RMSEs) of ice coverage were slightly decreased in Expt. 2 from the results in Expt. 1 during the winter of 2014–2015, except for Lake Ontario (Table 1). In the low-ice winter of 2015–2016, the RMSEs were similar among the three experiments. As anticipated, the simulated snow and ice conditions in Expt. 3 were nearly identical to that in Expt. 2. Therefore, in the following discussion in this section, we will focus on the results from Expts. 1 and 2 only.


| Superior | Huron | Michigan | Erie | Ontario | ||
|---|---|---|---|---|---|---|
| December 2014 to May 2015 | Control (Expt. 1) | 14.1 | 9.7 | 5.4 | 6.6 | 8.0 |
| Precipitation (Exp. 2) | 12.5 | 9.0 | 5.4 | 6.1 | 7.9 | |
| December 2015 to May 2016 | Control (Expt. 1) | 5.3 | 3.4 | 1.6 | 6.0 | 2.4 |
| Precipitation (Exp. 2) | 5.4 | 3.5 | 1.5 | 5.9 | 2.4 | |
Figure 5 shows the time series of simulated ice volume for each of the Great Lakes. A notable reduction in ice volume occurred by inclusion of precipitation (i.e., snow). The maximum ice thickness, which was obtained by dividing the ice volumes by the ice areas, was reduced by including precipitation (i.e., Expt. 2 versus Expt. 1) in all lakes for 2014–2015: 7, 6, 2, 8, and 1 cm for Lakes Superior, Huron, Michigan, Erie, and Ontario, respectively. In the low-ice winter of 2015–2016, the reduction in ice volume was less evident, but the maximum ice thickness was reduced by 3 and 1 cm in Lake Huron and Ontario, respectively. In the precipitation experiments, the shorter ice durations by 1–4 days were seen for most of the lakes (Table 2).

| Superior | Michigan | Huron | Erie | Ontario | ||
|---|---|---|---|---|---|---|
| 2014–2015 season | Control (Expt. 1) | 9 Jan. to 1 May | 8 Jan. to 11 Apr. | 6 Jan. to 27 Apr. | 8 Jan. to 20 Apr. | 18 Jan. to 11 Apr. |
| Precipitation (Exp. 2) | 9 Jan. to 30 Apr. | 9 Jan. to 10 Apr. | 6 Jan. to 23 Apr. | 8 Jan. to 17 Apr. | 18 Jan. to 9 Apr. | |
| NIC | 11 Jan. to 1 May | 5 Jan. to 1 May | 4 Jan. to 26 Apr. | 8 Jan. to 20 Apr. | 14 Jan. to 11 Apr. | |
| 2015–2016 season | Control (Expt. 1) | 14 Feb. to 12 Mar | 21 Jan. to 9 Mar. | 17 Jan. to 15 Mar. | 18 Jan. to 1 Mar. | — |
| Precipitation (Exp. 2) | 13 Feb. to 13 Mar | 21 Jan. to 9 Mar. | 17 Jan. to 15 Mar. | 18 Jan. to 1 Mar. | — | |
| NIC | 14 Feb. to 8 Mar. | 20 Jan. to 8 Mar. | 17 Jan. to 17 Mar. | 19 Jan. to 28 Feb. | 15 Feb. to 7 Mar. | |
- Note. Ice duration is defined as the period from when 5-day running mean ice coverage reaches 10% for the first time to when it goes below 10% for the last time in the season (i.e., December–May).
The modeled snow cover on the ice in the precipitation experiment (Expt. 2) was 0–8 cm in the winter of 2014–2015 and less than 1 cm in the winter of 2015–2016 (Figure 6). The observations of snow cover on the ice that can be compared with the model results are rare. There have been limited visual observations of snow depth on the ice in Lake Erie reported by the U.S. Coast Guard during their helicopter flyovers, in collaboration with NOAA Great Lakes Environmental Research Laboratory. These observations occurred on 1 or 2 days per year in late February or early March from 2008 to 2015 (with a few skipped years). In the winter of 2014–2015, the reported values on 25 and 26 February 2015 over five stations in Lake Erie ranged from 0 to 10 cm. The range does not conflict with the modeled values (Figure 6), but care should be taken as the range from the reports is for point values, while Figure 6 shows the lake-wide means. The ice surface albedo ranged from 6% (bare, very thin ice) to 60% in the control experiment. For comparisons, Bolsenga (1969) reported 10–46% for snow-free ice in Lake Superior, being consistent with the model results. The modeled albedo increased by 0–20% when precipitation was considered. This is consistent with Bolsenga (1987) who reported 62–80% for snow albedo in Michigan locations. By late March, snow melted away and the ice surface albedo went down quickly to be similar to that in Expt. 1 (Figure 6). The spatial patterns of snow depth presented high spatial heterogeneity (Figure 7), which also contributed to the spatial variation in the ice surface albedo.


3.2 Snow Impacts on Ice Conditions
Inclusion of snow reduced the total volume of ice despite that snow cover on the ice increased surface albedo in early winter of 2014–2015 (Figure 6), which reflected more sunlight back to the atmosphere and therefore prevented ice surface from melting. Note that snow also served as an input to the ice mass by formation of snow-ice. Apparently, the thermal insulation of snow cover on the ice (Sturm et al., 2002; Warren et al., 1999) overcame the impacts of increased surface albedo and snow-ice formation. In the precipitation experiment (Expt. 2), the surface albedo quickly reduced in March, when the atmosphere started to warm, and snow cover started to melt away. In this period, the surface albedo in the precipitation experiment was slightly below that in the control experiment (e.g., Lake Erie in the winter of 2014–2015, Figure 6). This was because the modeled lake ice was thinner in Expt. 2 than in Expt. 1. As described in section 2.1, the parametrized ice albedo reduced as ice became thinner (Briegleb, 1992). By including snow, the modeled ice thickness became thinner due to the dominating thermal insulation effect as described above. This resulted in the lower surface albedo once the ice lost snow cover on it. The shortened ice duration by inclusion of precipitation was likely because the reduced total ice volume by snow cover on the ice (manifestation of precipitation) allowed earlier delay of lake ice. Overall, inclusion of precipitation improved the models in reproducing ice areas and durations, as demonstrated by the reduction in RMSE of ice coverage.
In the version of model used in the study, the amount of formed snow ice was not directly outputted, not allowing the direct assessment. However, estimations of maximum possible lake-wide mean snow depth were made using the modeled freezing precipitation amount, that is, P in equation 2. The estimates resulted in 6, 3, 3, 2, 5 cm (4, 2, 5, 1, 2 cm) for Lake Superior, Lake Huron, Lake Michigan, Lake Erie, and Lake Ontario for the winter of 2014–2015 (2015–2016). If they were all to form snow ice, the estimated maximum possible additions to the total ice volume would be 1.7, 0.7, 0.7, 0.2, 0.3 km3 (1.2, 0.4, 1.1, 0.1, 0.1 km3). These numbers were fairly small compared with the range in Figure 5. Thus, inputs from snow ice formation were likely minor to the total ice volume. This is somewhat contradictory to the findings in the polar oceans (e.g., Fichefet & Morales Maqueda, 1999) and other lakes (e.g., Ohata et al., 2017), but in these study areas, snowfall amounts over the ice surface were much larger than the estimates in this study, providing favorable conditions for snow ice formation.
For comparison with classical freezing degree day (FDD) ice growth models, Tables 3 and 4 show the estimated maximum ice thickness by the equations proposed by Lebedev (1938), referred in Bilello, 1961) and Ashton (1989). Note that the Lebedev model was based on the sea ice observations in the Arctic shelf seas. The Ashton model was for lake ice but no snow cover was considered. Overall, in the winter of 2014–2015, the estimates by the simple FDD ice growth models were in agreement with the three-dimensional ice-hydrodynamic model results. The FDD ice growth models tended to provide higher ice thickness. This was pronounced in Lake Superior, Lake Huron, and Lake Michigan in the winter of 2015–2016 and in all lakes in the winter of 2016–2017. One obvious reason was heat (i.e., warm water) carried from the preceding summer, which prevented ice from growing in reality. An FDD ice growth model typically takes air temperature as the single input for its simplicity and therefore is limited in application to dimictic lakes like the Great Lakes. Snow inclusion in the Lebedev model resulted in some reductions of ice thickness estimates for a constant snow depth hs of 6 cm, as was seen in the three-dimensional model results. However, with hs = 3 cm, ice thickness estimates resulted in unrealistic increase from the no snow results in both winters. Again, the Lebenov model was an empirical model and the data fitted to the equation had hs and FDD mostly greater than 10 cm and 1,000°C/day, respectively (Figure 11 in Bilello, 1961). Therefore, its application to the Great Lakes should be limited as hs and FDDs are typically below these ranges in a normal winter. The Ashton model provided closer results to the three-dimensional model results than the Lebenov model did, but the higher ice thickness estimates were also pronounced in the winter of 2015–2016, likely due to not accounting of warm water carried from the preceding summer. Thus, the simple FDD ice growth models were useful for estimating ice thickness at the first order but tended to provide higher estimates than the three-dimensional model did due to lack of the thermal history in the lakes. It should also be noted that the three-dimensional model also takes account of mechanical thickening of ice, such as ridging and rafting.
| Superior | Huron | Michigan | Erie | Ontario | ||
|---|---|---|---|---|---|---|
| 2014–2015 | Maximum FDD (°C/day) | 1,008.5 | 809.2 | 533.5 | 618.0 | 582.2 |
| Lebedev (1938) | ||||||
| no snow, hi = 1.33FDD0.58 | 73.4 | 64.6 | 50.8 | 55.3 | 53.4 | |
| hs = 3 cm, hi = 1.45FDD0.62hs−0.15 | 76.9 | 67.1 | 51.8 | 56.8 | 54.7 | |
| hs = 6 cm, hi = 1.45FDD0.62hs−0.15 | 69.3 | 60.5 | 46.7 | 51.2 | 49.3 | |
| Asthon (1989) | ||||||
| no snow hi = [2κ/ρiLiFDD + (κ/Fs0)2]0.5 − κ/Fs0 | 68.7 | 60.6 | 47.5 | 51.8 | 50.0 | |
| 3D FVCOM-CICE model | ||||||
| no snow | 56 | 50 | 42 | 57 | 56 | |
| with snow | 49 | 44 | 40 | 49 | 55 | |
- Note. Maximum FDD for each winter was obtained as a cumulative value since the first day when daily air temperature went below 0°C. hs is a constant snow depth in cm. In the Ashton (1989) equation, FDD's unit should be converted to °C/s. κ is the thermal conductivity of ice (W/m2) and Fs0 is a constant value for the representative turbulent sensible heat flux (Fs0 = 20 W/m2 is used).
| 2015–2016 | Maximum FDD (°C/day) | 368.2 | 263.5 | 167.8 | 110.4 | 122.5 |
| Lebedev (1938) | ||||||
| no snow, hi = 1.33FDD0.58 | 40.9 | 33.7 | 26.0 | 20.4 | 21.6 | |
| hs = 3 cm, hi = 1.45FDD0.62hs−0.15 | 41.2 | 33.5 | 25.3 | 19.5 | 20.8 | |
| hs = 6 cm, hi = 1.45FDD0.62hs−0.15 | 37.1 | 30.2 | 22.8 | 17.6 | 18.8 | |
| Ashton (1989) | ||||||
| no snow hi = [2κ/ρiLiFDD + (κ/Fs0)2]0.5 − κ/Fs0. | 38.0 | 30.8 | 23.0 | 17.3 | 18.6 | |
| 3D FVCOM-CICE model | ||||||
| no snow | 11 | 24 | 21 | 11 | 15 | |
| with snow | 11 | 21 | 21 | 11 | 14 | |
3.3 Precipitation Impacts on Heat Fluxes and Water Temperature
The lake-wide mean precipitation heat fluxes Hlp and Hsp were in general intermittent due to their association with snow or rain events but presented expected seasonality (Figure 8). Snowfall during winter caused negative Hlp, (i.e., cooling of water surface). The daily values of Hlp were below −50 W/m2 even during major snowstorm events. This range was smaller than the net heat flux in the similar timeframe (−500–0 W/m2, Figure 9), however, not negligible. With a simple column equation, the water surface with a mixed layer depth of 10 m would be cooled at 0.1°C/day under a constant cooling of −50 W/m2. Depending on the winter mixed layer depth, which can be shallower (deeper) than 10 m at low (high) wind conditions, the rate could be faster or slower. The sensible heat flux due to precipitation, Hsp, was overall smaller in terms of magnitude than Hlp. Positive Hsp was due to warm rain relative to the water surface from spring to summer.


Water surface temperature was reasonably simulated by the models in comparison with the analysis from GLSEA (Figure 10a). The contrast between the anomalously cold and warm winters was demonstrated in the minimum values of lake-wide mean water surface temperature: In the winter of 2014–2015, the minimum values were around the freezing temperature for all five lakes, while in the winter of 2015–2016, the minimum daily values ranged from 0.3°C (Lake Erie) to 2.3°C (Lake Ontario). Inclusion of precipitation resulted in only minor changes in lake-wide mean water surface temperature. Biases in the modeled lake-wide mean water surface temperature were very similar among the three experiments (Figure 10b). However, the three experiments presented minor differences in three occasions. First, although not clearly visible in the lake-wide mean time series in Figure 10b, major snowstorms caused intermittent, local cooling at the water surface. Figures 11a and 11b show an example of such cooling on 3 February 2015, after a major snowstorm over the Great Lakes region resulted in significant Hlp particularly over Lake Ontario and Lake Michigan. However, this signal was buried in the lake-wide average time series (Figure 10). Second, from late April to early May in 2015, there was slight warming in the two precipitation experiments relative to the control experiment (Expt. 1) in Lake Superior, Lake Huron, and Lake Erie. This coincided with the increase in the incoming shortwave radiation and earlier decay of ice cover in the precipitation experiments during the same time (Figure 9b and Table 2). The warming in Expts. 2 and 3 during this period was likely because the 1–4 days earlier ice-off allowed shortwave radiation to reach the water surface and therefore warmed the water temperature earlier. The warming was most notable in Lake Erie, Lake Superior, and part of Lake Huron, as in example snapshots on 16 April 2015 (Figures 11c and 11d). This impact was not directly due to Hlp or Hsp, but the indirect consequence of snow cover on the ice leading to the earlier decay of the ice. Third, from February 2015 to May 2016, the water surface temperature was slightly cooled in the precipitation experiments. This weak signal can be seen in the reduced biases for Lake Superior and Huron (Figure 10b). This coincides with the negative Hlp during the period (Figure 8a). This feature was amplified in Expt. 3, where the response to the precipitation perturbation was constrained to water temperature and therefore maximized. The cooling was most evident in Lake Superior, Lake Huron and Lake Michigan, as in example snapshots on 30 May 2016 (Figures 11e and 11f). The lake-wide mean water surface temperature could have been more sensitive to Hlp in the winter of 2015–2016 due to the low ice coverage during the season. The amplitudes of cooling in Figure 11 appear to be at the same order as the column model estimate by Hlp provided earlier in this study. However, even so, the Hlp impacts on the water surface temperature appear to be small compared with the biases against the observation and the seasonal variability, given that these signals are hard to see in the lake-wide mean time series in Figure 10.


The RMSEs of water surface temperature were improved in the precipitation experiments, but only slightly (Table 5). Thus, the overall precipitation impacts on the water surface temperature were found to be local and intermittent and be very minor from a lake-wide average perspective.
| Superior | Huron | Michigan | Erie | Ontario | ||
|---|---|---|---|---|---|---|
| December 2014 to November 2015 | Control (Expt. 1) | 0.69 | 0.63 | 0.85 | 0.72 | 0.99 |
| Precipitation (Exp. 2) | 0.69 | 0.63 | 0.84 | 0.70 | 0.99 | |
| Precipitation (Exp. 3) | 0.69 | 0.63 | 0.83 | 0.69 | 0.99 | |
| December 2015 to November 2016 | Control (Expt. 1) | 1.05 | 0.77 | 0.79 | 0.58 | 1.09 |
| Precipitation (Exp. 2) | 1.05 | 0.63 | 0.78 | 0.58 | 1.08 | |
| Precipitation (Exp. 3) | 1.03 | 0.62 | 0.77 | 0.57 | 1.1 | |
3.4 Model Uncertainties
Major uncertainties in this study were from overlake precipitation estimates and representation of snow cover on the ice. Even though we used the “best” available estimates of overlake precipitation, we consider that the observations included higher uncertainty than the other meteorology variables due to its extremely high variability in space and time. The threshold to distinguish snow and rain (i.e., the wet-bulb temperature −2°C, see section 2.1) was empirical (Knox et al., 2017; Olsen, 2003) although it should serve as a first-order accurate threshold at least. It should be noted that the current HRRR product provides snowfall. However, the authors did not have access to the variable in the archive to cover the entire simulation period.
Like in many other models, representation of snow cover on the ice accompanied various approximations in this study. These include snow blown by winds, aging of snow, and the associated changes in the density and heat conductivity of snow. The snow-ice formation parameterization is based on the hydrostatic balance and highly simplified. In reality, the flooding of snow-ice interface can occur more frequently in a wave environment, which could allow more formation of snow ice. Relative importance of these processes is still unclear, not only in the Great Lakes but also in the polar regions (Lei et al., 2016, 2017; Webster et al., 2018). In order to determine which processes are significant in controlling Great Lakes ice cover (and therefore a priority for modeling), more mature observations are needed to examine models. In particular, spatial observations for ice thickness, snow depth on the ice, and snow/ice albedo by satellite or aircraft missions would be valuable to aid in understanding lake ice mass, energy budget, and the roles of snow in these components. Process-oriented observations would also greatly help to inform the development of lake ice models, including both numerical models and simple FDD ice growth models.
4 Summary and Conclusions
We evaluated the impacts of precipitation on the Great Lakes ice cover and water temperature using the state-of-the-art ice-hydrodynamic model, with regard to snow accumulation on the ice and the heat fluxes associated with rain and snow. The ice-hydrodynamic simulations were conducted for two anomalously cold and warm winters using the coupled FVCOM-CICE model for all of the Great Lakes. Modeled ice coverage and water surface temperature reasonably captured the seasonal and interannual variations in the observations. In both winters, the total ice volumes were reduced by the inclusion of snow cover in all five lakes, resulting in thinner lake ice and 1–4 days earlier ice-off. The earlier ice-off in the precipitation experiments resulted in earlier warming of water surface temperature, particularly in Lake Erie in the winter of 2014–2015, because the incoming shortwave radiation was allowed to reach the water surface earlier. Snow ice formation was unlikely a major input tot the total ice volume, largely due to the low overlake snowfall, but the relatively large uncertainty in overlake precipitation estimates posed some ambiguity on this conclusion. Estimates of maximum ice thickness were consistent among those from the coupled FVCOM-CICE model and the classical FDD ice growth models by Lebedev (1938) and Ashton (1989), but the classical FDD ice growth models were limited in representing impacts of snow impacts on ice growth and thermal history of the lakes.
Direct impacts of precipitation (snow or rain) on the water surface temperature were fairly minor in comparison with the seasonal variation of the water surface temperature. During the major snowstorm events, the open water areas received snow fall (i.e., a form of precipitation during late fall and winter) directly and cooled down by the release of latent heat due to snow melting (Hlp). The cooling due to Hlp was notable locally and intermittently during and after the snowstorm events, but its impacts on lake-wide mean temperature were not outstanding. The impacts of the precipitation sensible heat flux (Hsp) were negligible relative to Hlp or the other heat flux components (e.g., the turbulent heat fluxes).
The RMSEs for ice extent and lake-wide mean water surface temperature were slightly improved by the inclusion of precipitation, as a result of combined impacts described above (mostly from snowfall on the ice and water surface). The fact that snow cover inhibited ice growth and resulted in slightly earlier ice-off dates in the winter of 2014–2015 is compelling, reversing the commonly held notion that snow cover reflects the incoming shortwave radiation more than that of bare ice, which in turn extends the ice period. The findings from the model experiments during the two anomalously cold and warm winters indicate the combined impacts of precipitation described above may be expected in other normal years with average overlake air temperature and ice condition, but these impacts could be amplified or diminished depending on how much snow falls over the lakes in the winter of interest.
This study presented that winter precipitation, particularly snowfall, is an important factor in the winter energy budget over ice and water in the Great Lakes. In the changing climate, every component of the system is expected to change in a complexed way: Future projections suggest more precipitation on average, possibly leading to increased snow fall over the lakes, but rising temperatures will cause more winter precipitation to fall as rain as opposed to snow across the region by late century (Notaro, Bennington, & Lofgren, 2015; Notaro, Bennington, & Vavrus, 2015). Great Lakes ice cover is expected to continue to decrease, but there remains strong year-to-year variability, and high ice years are still possible (Mason et al., 2016; Wang et al., 2012). More observational studies are needed for overlake precipitation, snow cover, albedo, and ice thickness to reduce model uncertainties and ultimately improve the regional climate projections.
Acknowledgments
The study did not generate any new observational data. The merged CaPA-MPE precipitation data set can be accessed at the website of the Midwest Regional Climate Center (https://mrcc.illinois.edu/gismaps/naprecip.htm). The analyses of the Great Lakes ice coverage and surface water temperature can be accessed at NOAA CoastWatch Great Lakes node website (https://coastwatch.glerl.noaa.gov/statistic/statistic.html). The HRRR surface meteorological data set can be accessible at the University of Utah HRRR Download page (http://home.chpc.utah.edu/~u0553130/Brian_Blaylock/cgi-bin/hrrr_download.cgi) and NOAA National Centers for Environmental Prediction (https://www.nco.ncep.noaa.gov/pmb/products/hrrr/). The original source codes for the Finite Volume Community Ocean Model (FVCOM) can be accessible at http://fvcom.smast.umassd.edu/FVCOM/Source/code.htm (user registration required). The updated source codes of FVCOM used in this study, the namelist files, and the representative model results are archived at the Deep Blue repository at the University of Michigan with the data set title “FVCOM-UGCICE source codes, namelists, and model results for assessment of precipitation impacts on Great Lakes ice cover and water temperature” (https://doi.org/10.7302/4mpr-4364). This research was carried out with support of the National Oceanic and Atmospheric Administration (NOAA) Great Lakes Environmental Research Laboratory (GLERL), awarded to the Cooperative Institute for Great Lakes Research (CIGLR) through the NOAA Cooperative Agreement with the University of Michigan (NA12OAR4320071). This is GLERL contribution 1947 and CIGLR contribution 1161. We thank Yin Min Goh at the University of Michigan for preliminary data analysis for the precipitation heat flux and Gregory Lang for obtaining the atmospheric forcing data from the HRRR archive.